In this study, we have steroid information of 31 male patients at three different time points after trauma (<1h, 4-12h, 24-48h). Only 26 have complete steroid information for the three time points. We had also immunological information for these patients but it was taken out from the study. Clinical information such as age, MODS, LOS, ISS, Drugs, Fluids.., is also included. 34 healthy controls are also included for reference.

Two analysis with different steroid’s concentration have been performed given we had two different LOQ values given. The first measurements taken with old LOQ and the new revalidated LOQ.

This is the one with the REVALIDATED LOQ

#Loading packages
#install_github("kassambara/easyGgplot2")
require(readxl)
require(reshape2)
require(plyr)
require(lubridate)
require(stringr)
require(ggrepel)
require(data.table)
require(ggpubr)
require(readr)
require(tidyverse)
require(dplyr)
require(knitr)
require(htmltools)
require(factoextra)
require(Metrics)
require(e1071)
require(ggplot2)
require(gridExtra)
require(brglm)
require(ggstatsplot)
require(compareGroups)
require(knitr)
require(kableExtra)
require(gtsummary)
require(dunn.test)
require(broom)
require(reshape2)
require("patchwork")
require(gam)
require(mgcv)
require(ComplexHeatmap)

#this might need github too
#packageurl <- "https://cran.r-project.org/src/contrib/Archive/repr/repr_1.0.2.tar.gz"
#install.packages(packageurl, repos=NULL, type="source")

set.seed(132)

1 Data import

1.1 Reference steroid’s concentration (LOQ)

File with LOQ concentrations (lowest concentration at which the analyte can not only be reliably detected).It includes the parameter we used at the start, and at which the original analysis was developed (original) and the corrected/revalidated one, where many datapoints were filtered as a result.

Again this file contains all analysis performed with REVALIDATED LOQ

Steroids <- read_xlsx("SteroidsLow.xlsx")
Steroids <- Steroids[complete.cases(Steroids),c(2,3,4)]
names(Steroids) <- c("variable", "LOQ_Original", "LOQ_Revalidation")

1.2 Healthy Patients (controls)

34 healthy controls with steroid information. Standard ranges were also recovered and measurements outside those ranges were taken out.

#outliers per steroid. 
#saveRDS(Names, "names.rds")

RealHealthy <- function(HealthyLong){
  
Names <- data.frame(variable = unique(filter(HealthyLong, Descriptor == "Steroid")$variable))
Names["LimitLow"] <- c("30","10","1","0.7","2","0.1","2","0.6","1","0.09","0.3","2","0","0.07","0","1", rep("0",7)) #said 1 in progesterone (second to last but then all readings would be off)
Names["LimitHigh"] <- c("100","650","2.6","3.9","13","0.7","40","4","7","0.6","40","35","10","2.4","10","20",rep("100",7)) #much less in170HP but then off too * also 11KA4 very high numbers in general - muh than trauma and reference
Names

a <- full_join(HealthyLong,Names)

a$LimitLow <- as.numeric(a$LimitLow)
a$LimitHigh <- as.numeric(a$LimitHigh)

a1 <- a %>% 
  mutate(TakeOut = value >= (LimitLow)  & value <= (LimitHigh)) #take out ten percent of readings? -- Yes

HealthyAll <- a1 %>% 
  filter( TakeOut == "TRUE") %>% 
  select(-c("LimitLow","LimitHigh", "TakeOut"))

return(HealthyAll)

}

Healthy <- read_delim("FinalHealthyPat.csv", 
    ";", escape_double = FALSE, locale = locale(decimal_mark = ",",
 grouping_mark = "."), trim_ws = TRUE) %>% 
  select(-Age) #Take out Age for healthy - should be age matched either way.

names(Healthy)<- c("ID", "Label", "Cortisone", "Cortisol", "11-KTestosterone","11-KA4", "11β-OHA4","11β-OHT", "Corticosterone", "11-deoxycortisol", "Androstenedione", "11-deoxycorticosterone", "Testosterone", "DHEA", "17-hydroxyprogesterone", "Dihydrotestosterone", "Progesterone", "DHEAS", "Cortisol/Cortisone", "DHEA/DHEAS", "Cortisol/DHEAS", "Cortisol/DHEA", "DHEA/Testosterone", "Testosterone/Cortisol", "Aldosterone") #distinguish between trauma and healthy control

HealthyLong <- reshape2::melt(Healthy, id.vars=c("ID","Label"))

HealthyLong <- HealthyLong %>% 
  mutate(Time = "HC", 
         Descriptor = "Steroid") 

HealthyLong <- RealHealthy(HealthyLong) #LOQ feeds out of this one so all fine

HealthyLongLOQ <- full_join(HealthyLong,Steroids) %>%
  replace_na(list(LOQ_Original = 0, LOQ_Revalidation = 0)) %>% 
  filter(value >= LOQ_Revalidation)

HealthyWideLOQ <- HealthyLongLOQ %>% 
  select(ID, Label, variable, value, Descriptor) %>%
  pivot_wider(names_from = variable, values_from = value)

#take out healthy patients that do not follow normal ranges (??)

1.3 Non healthy Patients (trauma)

31 trauma patients with cytokine, immunological and clincal information at three different time points.

(Lesson: pivot longer will put problems with types and classes of columns but melt wont).

Trauma <- read_csv("FinalDatasetAprilwCytokines.csv") #non healthy patients

names(Trauma)[1] <- "ID"

TraumaLong <- reshape2::melt(Trauma, id.vars=c(1))

TakeOut <- c(grep("^[*]",TraumaLong$value),which(TraumaLong$value %in% c("ND","OOR <"))) #Take out spurious values 

TraumaLong <- TraumaLong[-c(TakeOut),] %>% 
  drop_na() %>% 
  mutate(Label = "Trauma")

TraumaLong[,c("variable","Time") ] <- str_split_fixed(TraumaLong$variable, "-", 2)

TraumaLong$Time <- gsub(" ","",TraumaLong$Time)

TraumaLong$Time <- as.factor(TraumaLong$Time)

levels(TraumaLong$Time) <- c("Clinical","T1","T2","T3")

'%ni%' <- Negate('%in%')

TraumaLong <- TraumaLong %>%
  filter(variable %ni% c("IL1ra","IL6", "IL8","IL10","GCSF", "MCP1", "TNFa","IL12","IL17", "Lymphocyte (106/L)","Monocyte (106/L)","Neutro:Lymph ratio","Neutrophil (106/L)","WBC (106/L)","IG (106/L)")) %>%
  mutate(Descriptor = ifelse(Time == "Clinical", "Clinical", "Steroid"))

TraumaLong$variable <- as.factor(TraumaLong$variable)

levels(TraumaLong$variable) <-c("11-deoxycortisol","11-KA4","11-KTestosterone","11β-OHA4","11β-OHT","17-hydroxyprogesterone", "Androstenedione","Age","Alcohol", "Aldosterone","Anaesthetic","Corticosterone","Cortisol",
"Cortisol/Cortisone","Cortisol/DHEA","Cortisol/DHEAS","Cortisone","DHEA","DHEA/DHEAS","DHEA/Testosterone","DHEAS",
                              "Dihydrotestosterone","11-deoxycorticosterone","Fluids","GCS","GoldenHour","ISS","MODS", "Lactate","LOS","Mechanisms","NISS","Outcome","Drug","Progesterone","Sex","Testosterone","Testosterone/Cortisol",
                                 "T1","T2","T3")


TraumaWide <- TraumaLong %>%
  select(- c(Label, Descriptor)) %>%
  pivot_wider(values_from = value, names_from= c(Time,variable))


TraumaLongLOQ <- full_join(TraumaLong,Steroids) %>%
  replace_na(list(LOQ_Original = 0, LOQ_Revalidation = 0)) %>% 
  filter(value >= LOQ_Revalidation)

names(TraumaLongLOQ)[1] <- "ID"

TraumaWideLOQ <- TraumaLongLOQ %>% 
  select(ID, Label, Time, variable, value) %>%
  pivot_wider(names_from = c(Time,variable), values_from = value)

2 Data summary (Revalidated LOQ)

CreateTablesTP <- function(TraumaWide,TraumaWideLOQ, TimePointOut1,TimePointOut2, TimePointIn ) {
  
  TraumaT1 <- TraumaWide %>% 
  select(-c(contains(TimePointOut1), contains(TimePointOut2), "ID")) %>% 
  #drop_na() %>%
  mutate_at(vars(starts_with(TimePointIn), "Clinical_Age", "Clinical_GCS", "Clinical_ISS", "Clinical_NISS", "Clinical_LOS", "Clinical_Lactate", "Clinical_GoldenHour" ),funs(as.numeric))

names(TraumaT1) <- str_remove_all(names(TraumaT1), TimePointIn)

TableTraumaT1 <- TraumaT1 %>%
  tbl_summary( 
              type = list(all_numeric() ~ "continuous",
                          all_logical() ~ "categorical") ,
              missing_text = "(Missing)"
  ) %>%
  bold_labels() %>%
  modify_header(label = "**Variable**") 

#

TraumaT1LOQ <- TraumaWideLOQ %>% 
  select(-c(contains(TimePointOut1), contains(TimePointOut2), "ID")) %>% 
  #drop_na() %>%
  mutate_at(vars(starts_with(TimePointIn), "Clinical_Age", "Clinical_GCS", "Clinical_ISS", "Clinical_NISS", "Clinical_LOS", "Clinical_Lactate", "Clinical_GoldenHour" ),funs(as.numeric))

names(TraumaT1LOQ) <- str_remove_all(names(TraumaT1LOQ), TimePointIn)

TableTraumaT1LOQ <- TraumaT1LOQ %>%
  tbl_summary( 
              type = list(all_numeric() ~ "continuous",
                          all_logical() ~ "categorical") ,
              missing_text = "(Missing)"
  ) %>%
  bold_labels() %>%
  modify_header(label = "**Variable**") 

return(list(TableTraumaT1,TableTraumaT1LOQ ))
  
}


TP1 <- CreateTablesTP(TraumaWide,TraumaWideLOQ, "T2_","T3_", "T1_" ) 
TP2 <- CreateTablesTP(TraumaWide,TraumaWideLOQ, "T1_","T3_", "T2_" ) 
TP3 <- CreateTablesTP(TraumaWide,TraumaWideLOQ, "T2_","T1_", "T3_" ) 
  

##########################


HealthyTable <- Healthy %>% 
  #drop_na() %>%
  select(-c(ID, Label)) %>%
  tbl_summary( 
             type = list(all_numeric() ~ "continuous",
                        all_logical() ~ "categorical") ,
              missing_text = "(Missing)"
  ) %>%
  bold_labels() %>%
  modify_header(label = "**Variable**") 

HealthyTableLOQ <- HealthyWideLOQ %>% 
  #drop_na() %>%
  select(-c(ID, Label, Descriptor)) %>%
  tbl_summary( 
              type = list(all_numeric() ~ "continuous",
                          all_logical() ~ "categorical") ,
              missing_text = "(Missing)"
  ) %>%
  bold_labels() %>%
  modify_header(label = "**Variable**")

2.1 Table 1.

tbl_merge(tbls = list(HealthyTableLOQ, TP1[[2]], TP2[[2]], TP3[[2]]), tab_spanner = c("**Healthy Controls**", "**Time Point <1h**", "**Time Point 4-12h**", "**Time Point 24-48h**") )
Variable Healthy Controls Time Point <1h Time Point 4-12h Time Point 24-48h
N = 341 N = 312 N = 312 N = 312
Cortisone 62 (53, 70) 53 (42, 60) 57 (50, 70) 34 (22, 45)
(Missing) 3 2 5
Cortisol 273 (210, 335) 458 (386, 573) 510 (325, 605) 214 (114, 321)
(Missing) 2 5
11-KTestosterone 1.50 (1.16, 1.60) 1.4 (0.8, 2.4) 10 (1, 12) 0.8 (0.7, 11.8)
(Missing) 16 2 5
11-KA4 3.20 (2.65, 3.45) 1.8 (1.4, 2.8) 1.69 (1.12, 3.01) 1.50 (1.10, 3.97)
(Missing) 23 2 5
11β-OHA4 8.45 (6.32, 10.07) 16 (12, 22) 12 (8, 21) 9.6 (5.2, 12.4)
(Missing) 12 2 8
11β-OHT 0.600 (0.600, 0.625) 0.92 (0.78, 1.23) 0.79 (0.78, 1.61) 0.6886 (0.6886, 0.6886)
(Missing) 27 9 25 30
Corticosterone 7 (5, 15) 61 (49, 88) 39 (15, 76) 7 (3, 15)
(Missing) 2 8
11-deoxycortisol 1.10 (0.80, 1.68) 4.5 (3.0, 7.6) 5.5 (3.0, 7.7) 1.17 (0.66, 2.58)
(Missing) 8 1 6 15
Androstenedione 2.97 (2.62, 4.32) 3.98 (3.24, 5.69) 3.20 (2.15, 5.45) 3.12 (2.09, 4.13)
(Missing) 3 12
11-deoxycorticosterone 0.30 (0.30, 0.35) 0.74 (0.40, 1.64) 0.95 (0.81, 2.32) 0.77 (0.56, 0.98)
(Missing) 27 9 13 29
Testosterone 14.5 (11.5, 17.4) 11.9 (8.7, 13.9) 3.83 (2.63, 5.75) 1.7 (0.9, 4.5)
(Missing) 2 7
DHEA 19 (13, 21) 32 (24, 54) 18 (14, 33) 8 (3, 15)
(Missing) 5 2 5
17-hydroxyprogesterone 3.25 (1.83, 4.07) 5 (4, 7) 4 (2, 11) 1.37 (0.93, 2.69)
(Missing) 3 15
Dihydrotestosterone 1.60 (1.15, 1.87) 1.08 (0.65, 1.29) 0.60 (0.60, 0.77) 0.60 (0.60, 0.60)
(Missing) 8 2 6
Progesterone 0.4000 (0.4000, 0.4000) 0.76 (0.60, 1.26) 1.04 (0.80, 1.93) 0.344 (0.339, 0.346)
(Missing) 31 3 12 28
DHEAS 9.3 (6.4, 11.1) 13 (8, 20) 8 (6, 18) 14 (10, 18)
(Missing) 2 5
Cortisol/Cortisone 4.53 (3.88, 6.16) 8.87 (8.03, 10.14) 9.3 (5.8, 11.5) 5.58 (4.57, 7.35)
(Missing) 2 5
DHEA/DHEAS 2.24 (1.56, 3.74) 3.00 (1.73, 3.87) 2.04 (1.33, 3.54) 0.77 (0.22, 1.55)
(Missing) 2 5
Cortisol/DHEAS 30 (22, 44) 37 (26, 55) 55 (26, 86) 15 (8, 25)
(Missing) 1 2 5
Cortisol/DHEA 13.6 (10.4, 17.1) 13 (10, 19) 110 (56, 230) 22 (13, 43)
(Missing) 2 5
DHEA/Testosterone 1.34 (0.92, 2.02) 3.19 (2.03, 4.91) 4.2 (2.8, 10.6) 3 (2, 9)
(Missing) 2 5
Testosterone/Cortisol 0.049 (0.039, 0.065) 0.024 (0.019, 0.033) 0.009 (0.004, 0.018) 0.010 (0.004, 0.027)
(Missing) 2 5
Aldosterone 0.60 (0.60, 0.60) 0.74 (0.60, 1.07) 0.70 (0.60, 1.05) 0.60 (0.60, 0.60)
(Missing) 2 5
Label
Trauma 31 (100%) 31 (100%) 31 (100%)
Clinical_Age 26 (22, 30) 26 (22, 30) 26 (22, 30)
Clinical_Sex
Male 31 (100%) 31 (100%) 31 (100%)
Clinical_GCS 15.00 (14.00, 15.00) 15.00 (14.00, 15.00) 15.00 (14.00, 15.00)
Clinical_ISS 16 (10, 21) 16 (10, 21) 16 (10, 21)
Clinical_NISS 22 (16, 29) 22 (16, 29) 22 (16, 29)
Clinical_Fluids 9 (29%) 9 (29%) 9 (29%)
Clinical_Anaesthetic 12 (39%) 12 (39%) 12 (39%)
Clinical_Drug
Cocaine 1 (3.2%) 1 (3.2%) 1 (3.2%)
Fentynl 1 (3.2%) 1 (3.2%) 1 (3.2%)
Ketamine 6 (19%) 6 (19%) 6 (19%)
Morphine 10 (32%) 10 (32%) 10 (32%)
NO 11 (35%) 11 (35%) 11 (35%)
Paracetemol 1 (3.2%) 1 (3.2%) 1 (3.2%)
Rocuromium 1 (3.2%) 1 (3.2%) 1 (3.2%)
Clinical_MODS 13 (42%) 13 (42%) 13 (42%)
Clinical_Mechanisms
Assault 1 (3.2%) 1 (3.2%) 1 (3.2%)
Assault and Penetrating 4 (13%) 4 (13%) 4 (13%)
Blunt Impact 2 (6.5%) 2 (6.5%) 2 (6.5%)
Fall From Height 1 (3.2%) 1 (3.2%) 1 (3.2%)
Penetrating 8 (26%) 8 (26%) 8 (26%)
RTC 15 (48%) 15 (48%) 15 (48%)
Clinical_Alcohol 5 (16%) 5 (16%) 5 (16%)
Clinical_Outcome
Alive 31 (100%) 31 (100%) 31 (100%)
Clinical_LOS 10 (6, 20) 10 (6, 20) 10 (6, 20)
Clinical_Lactate 3.20 (2.60, 4.78) 3.20 (2.60, 4.78) 3.20 (2.60, 4.78)
Clinical_GoldenHour 36 (26, 46) 36 (26, 46) 36 (26, 46)
Clinical_T1
Afternoon 13 (42%) 13 (42%) 13 (42%)
Evening 7 (23%) 7 (23%) 7 (23%)
Morning 4 (13%) 4 (13%) 4 (13%)
Night 7 (23%) 7 (23%) 7 (23%)
Clinical_T2
Afternoon 12 (41%) 12 (41%) 12 (41%)
Evening 4 (14%) 4 (14%) 4 (14%)
Morning 12 (41%) 12 (41%) 12 (41%)
Night 1 (3.4%) 1 (3.4%) 1 (3.4%)
(Missing) 2 2 2
Clinical_T3
Afternoon 2 (7.7%) 2 (7.7%) 2 (7.7%)
Morning 24 (92%) 24 (92%) 24 (92%)
(Missing) 5 5 5

1 Median (IQR)

2 n (%); Median (IQR)

We have more missing values in revalidated than original, because this new revalidation has caused many measurements to become invalid (missing). To start with we had 5 missing patient values in T3 and 2 for T2.

AllDataLOQLong <- rbind(TraumaLongLOQ, HealthyLongLOQ)
AllDataLong <- rbind(TraumaLong, HealthyLong)

2.2 Boxplot steroid description

General steroid concentration changes at the three different time points

2.2.1 Steroid concentration figures

(k3[["Progesterone"]]   |  k3[["17-hydroxyprogesterone"]]) / 
  (k3[["11-deoxycortisol"]]|  k3[["Cortisol"]]) / 
 (k3[["Cortisone"]] | k3[["Cortisol/Cortisone"]])
Figure.1

Figure.1 

k3[["11-deoxycorticosterone"]] | k3[["Corticosterone"]] | k3[["Aldosterone"]]
Figure.2

Figure.2 

(k3[["DHEA"]] | k3[["DHEAS"]] | k3[["DHEA/DHEAS"]]) /
(k3[["Androstenedione"]] | k3[["Testosterone"]] | k3[["Dihydrotestosterone"]])
Figure.3

Figure.3 

(k3[["11β-OHA4"]]  | k3[["11β-OHT"]]) /
(k3[["11-KA4"]]  | k3[["11-KTestosterone"]] )
Figure.4

Figure.4

2.3 Golden Hour Study

Study of steroid concentration at the first time point (<1h). Use of Generalized Additive Models (GAMs) to extract trajectories through time.

2.3.1 Golden Hour GAMs

wrap_plots(pQuanta)   + plot_layout(guides = "collect")
Golden Hour GAM

Golden Hour GAM

Original plots with patients above/below normal healthy concentrations. GAMs tries to build two different models (around blue dots and orange dots). As you can see, those steroids without model are those where the separation between orange and blue is practically non existent (e.g DHEA/DHEAS)

wrap_plots(pQuant2a)  + plot_layout(guides = "collect")
Golden Hour All

Golden Hour All

2.4 Heatmaps

Mineralocorticoids <- c("11-deoxycorticosterone","Corticosterone","Aldosterone")
GlucocorticoidPrec <- c("Progesterone","17-hydroxyprogesterone","11-deoxycortisol")
Glucocorticoids <- c("Cortisol", "Cortisone" )
Androgens <- c("DHEA","DHEAS","Androstenedione")
OxygenatedAndro <- c("Testosterone", "Dihydrotestosterone","11-KA4","11-KTestosterone","11B-OHA4","11B-OHT")
Ratios <- c("Cortisol/Cortisone", "Cortisol/DHEA","Cortisol/DHEAS","DHEA/DHEAS","DHEA/Testosterone","Testosterone/Cortisol")

aa <- list(Mineralocorticoids,GlucocorticoidPrec, Glucocorticoids,Androgens, OxygenatedAndro, Ratios) 
names(aa) <- c("Mineralocorticoids","GlucocorticoidPrec", "Glucocorticoids","Androgens", "OxygenatedAndro", "Ratios")

ll <- data.frame(aa) %>% 
  pivot_longer(everything()) %>% 
  unique() %>%
  as.data.frame()

https://github.com/kassambara/ggpubr/issues/65 p value adjust

# it was working before and now not - issues with names. 

scale_this <- function(x){
  (x - mean(x, na.rm=TRUE)) / sd(x, na.rm=TRUE)
}


Heatmapping <- function(ChosenData,ll){

AverageHC <-  ChosenData %>% 
  filter(Time =="HC") %>%
  group_by(variable) %>%
  mutate(MeanValue = mean(as.numeric(value))) %>% 
  select(variable, MeanValue) %>%
  unique() 

#AverageHC$variable <- str_replace(AverageHC$variable,  "11β-OHA4", "11B-OHA4")
#AverageHC$variable <- str_replace(AverageHC$variable,  "11β-OHT", "11B-OHT")

MixingHM <- inner_join(ChosenData, AverageHC) #adds a column with mean HC value

FinalDataHM <- MixingHM %>% 
  filter( Descriptor %in% c("Steroid"), Time !="HC") 

FinalDataHM$variable <- str_replace(FinalDataHM$variable,  "11β-OHA4", "11B-OHA4")
FinalDataHM$variable <- str_replace(FinalDataHM$variable,  "11β-OHT", "11B-OHT")

FinalDataHM$variable <- droplevels(as.factor(FinalDataHM$variable))
FinalDataHM$value <- as.numeric(FinalDataHM$value)

FinalDataHM <- FinalDataHM %>%
  mutate(valueWHC = (value-MeanValue)) %>%
  group_by(variable) %>%
  mutate(value_scaled = scale_this(valueWHC)) %>%
  select(ID, variable, value_scaled, Time) %>%
  pivot_wider(id_cols = ID, names_from = c(variable, Time), values_from = value_scaled) 

Clinical <- ChosenData %>% 
  filter( Descriptor !=  c("Steroid"), Time !="HC") %>% 
  select(ID, variable, value, Time) %>%
  pivot_wider(id_cols = ID, names_from = c(variable, Time), values_from = value) %>% 
  as.data.frame()

FinalDataHMID <- FinalDataHM$ID
FinalDataHMMatrix <- as.matrix(FinalDataHM %>% 
                      select(-c(ID, '11B-OHT_T3','11B-OHT_T2'))) ## too many missing values
rownames(FinalDataHMMatrix) = FinalDataHMID

ll3 <- data.frame(str_split(colnames(FinalDataHMMatrix), "_", simplify = TRUE))
names(ll3) <- c("value", "time")
ll2 <- inner_join(ll3,  ll) 


annotation_col = data.frame(
   
    Time = ll2$time,
    Type  = ll2$name 
    #Mean = ll2a$Mean
)
rownames(annotation_col) = colnames(FinalDataHMMatrix)


annotation_row = data.frame(
    
    Age= as.numeric(Clinical$Age_Clinical),
    LOS = as.numeric(Clinical$LOS_Clinical),
    ISS = as.numeric(Clinical$ISS_Clinical),
    Fluids = as.factor(Clinical$Fluids_Clinical),
    Anaesthesia = as.factor(Clinical$Anaesthetic_Clinical)
)
rownames(annotation_row) = Clinical$ID


ann_colors = list(
  
  Type = c(Androgens = "light blue", GlucocorticoidPrec = "yellow", Glucocorticoids = "orange", Mineralocorticoids = "springgreen4",OxygenatedAndro= "dark blue", Ratios ="snow2"),
  Time = c(T1 = "green", T2 = "orange", T3 = "red"), 
  Anaesthesia = c(NO = "lawngreen", YES = "red2"), 
  Fluids = c(NO = "lawngreen", YES = "red2")#, 
  #Age =  c("white", "firebrick"),
  #ISS = c("white", "firebrick"),
  #LOS = c("white", "firebrick")
)


#cutree_cols  = 5,cutree_rows  = 3,
Heatmap <- pheatmap(FinalDataHMMatrix,  annotation_col = annotation_col, annotation_row = annotation_row, cutree_cols  = 5,cutree_rows  = 3,annotation_colors = ann_colors, fontsize_col = 9, labels_row =paste0("Patient.", 1:31))#, fontsize_row = 0


   # annotation_colors = ann_colors)
return(list(FinalDataHM, ll,Clinical, Heatmap ))
}

Once again, main difference with original LOQ data is the amount of missing values (grey cells) [ Had to delete ‘11B-OHT_T3’,‘11B-OHT_T2’ as had only one value left]

AldNo <- AllDataLOQLong %>% 
  filter(variable != "Aldosterone")
Heatmapping(AldNo,ll)[[4]]
Heatmap with NO LOQ change

Heatmap with NO LOQ change 

pdf("HeatmapAll.pdf", 17, 11)
Heatmapping(AldNo,ll)[[4]] #too many missing values!
dev.off()
## quartz_off_screen 
##                 2

2.4.1 General average

# it was working before and now not - issues with names. 

scale_this <- function(x){
  (x - mean(x, na.rm=TRUE)) / sd(x, na.rm=TRUE)
}



Average <-  AldNo %>% 
  filter(Time != "Clinical") %>%
  group_by(variable) %>%
  mutate(value_scaled = scale_this(as.numeric(value))) %>%
  ungroup() %>%
  group_by(variable,Time) %>%
  mutate(MeanValue = mean(as.numeric(value_scaled))) %>% 
  select(variable, MeanValue) %>%
  unique() 

Average$variable <- str_replace(Average$variable,  "11β-OHA4", "11B-OHA4")
Average$variable <- str_replace(Average$variable,  "11β-OHT", "11B-OHT")

AverageMat <- pivot_wider(Average,names_from = Time, values_from =MeanValue ) 

FinalDataHMID <- AverageMat$variable
AverageMatFinal <- AverageMat %>% as.data.frame() %>%
  dplyr::select(c("HC","T1", "T2", "T3"))
FinalDataHMMatrix <- as.matrix(AverageMatFinal)
rownames(FinalDataHMMatrix) <- FinalDataHMID

ll2 <- ll[match(rownames(FinalDataHMMatrix), ll$value),] %>% 
  filter(value != "Aldosterone") 

annotation_row = data.frame(
    
    Type  = ll2$name
)


rownames(annotation_row) = rownames(FinalDataHMMatrix) 


ann_colors = list(
  
  Type = c(Androgens = "light blue", GlucocorticoidPrec = "yellow", Glucocorticoids = "orange", Mineralocorticoids = "springgreen4",OxygenatedAndro= "dark blue", Ratios ="snow2")
)

######
AveragePlot <- pheatmap(FinalDataHMMatrix,  annotation_row = annotation_row, annotation_colors = ann_colors, fontsize_col = 9,cluster_cols=FALSE)#,

AveragePlot
Heatmap average

Heatmap average 

pdf("HeatmapAverageConce.pdf",  7,  5)
print(AveragePlot)
dev.off()
## quartz_off_screen 
##                 2